// P1Diffractor.cpp: implementation of the P1Diffractor class.
//
//////////////////////////////////////////////////////////////////////

#include "P1Diffractor.h"

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

P1Diffractor::P1Diffractor(char *DiffractorFile, 
						   Position *W1, Position *W2, Position *W3, char **atomLabel, int nAtom)
: Diffractor(DiffractorFile, W1, W2, W3, atomLabel, nAtom)
{
  newM(v1, NIntensity, NAtom);
  newM(v2, NIntensity, NAtom);
}


P1Diffractor::~P1Diffractor()
{
  deleteM(v1);
  deleteM(v2);
}


//////////////////////////////////////////////////////////////////////
// functions
//////////////////////////////////////////////////////////////////////

/// @ IDerivate input derivatives of chi-sq with respect to intensities
/// @ XDerivative output derivatives of chi-sq with respect to orthonormal atom coord.
///               
void P1Diffractor::TransformDerivatives(double *IDerivative, 
										Position *XDerivative)
{
  int    n,m;
  double x1, x2, s1, s2, s3, s4, s5, s6;

  for(n=1;n<=NAtom;n++)	
	  for(m=1;m<=3;m++)
		  FDeriv[n][m] = 0.0;


  double *w1,*w2;

  // calculate first derivatives with respect to fractional atom coord and
  // store these in FDeriv
  for(m=1;m<=NIntensity;m++)
	{
    // deriv of chi-squared with respect to Im * real part of structure factor 
	  x1 = -IDerivative[m]*ReAmplitude[m];
    // deriv of chi-squared with respect to Im * imaginary part of structure factor 
	  x2 =  IDerivative[m]*ImAmplitude[m];

    
	  s1 = DH[m] * x1; // H * x1
	  s2 = DH[m] * x2; // H * x2
	  s3 = DK[m] * x1; // K * x1
	  s4 = DK[m] * x2; // K * x2
	  s5 = DL[m] * x1; // L * x1
	  s6 = DL[m] * x2; // L * x2

	  w1 = v1[m]; // as a reminder v1[m] is number of atoms vector of numbers form factor * sin(2*PI(hx+ky+lz))
	  w2 = v2[m]; // form factor * cos(2*PI(hx+ky+lz))

    // in FDeriv store derivs of chi-sq with respect to fractional atom coord.
	  for(n=1;n<=NAtom;n++)
    {
		  FDeriv[n][1] += w1[n]*s1 + w2[n]*s2; 
		  FDeriv[n][2] += w1[n]*s3 + w2[n]*s4; 
		  FDeriv[n][3] += w1[n]*s5 + w2[n]*s6; 		
		}
	}

  // Convert FDeriv into XDerivative
  // which correspond to converting from fractional atom coor to orthonormal atom coor
  MakeXDerivative(XDerivative);	// multiply somewhere by 4*PI?
}


void P1Diffractor::FourierTransform(Position *AtomPosition)
{
  int    n,m;
  double SinKYplusLZ,CosKYplusLZ;

  double s,c;

  int h,k,l;

  for(n=1;n<=NIntensity;n++)
	{
	  ReAmplitude[n] = 0.0;	
	  ImAmplitude[n] = 0.0;
	}

  for(m=1;m<=NAtom;m++)
	{
	  SetUpTrig(AtomPosition[m]);

  	for(n=1;n<=NIntensity;n++)
		{
	  	h = H[n];
  		k = K[n];
  		l = L[n];
      
      // cal sin(2*PI(ky+lz)) = sin(2*PI(ky))*cos(2*PI(lz))+cos(2*PI(ky))*sin(2*PI(lz))
	  	SinKYplusLZ	= SinY[k]*CosZ[l] + CosY[k]*SinZ[l];
      // cal cos(2*PI(ky+lz)) = cos(2*PI(ky))*cos(2*PI(lz))-sin(2*PI(ky))*sin(2*PI(lz))
		  CosKYplusLZ = CosY[k]*CosZ[l] - SinY[k]*SinZ[l];

      // cal sin(2*PI(hx+ky+lz))
		  s = SinX[h]*CosKYplusLZ + CosX[h]*SinKYplusLZ;
      // cal cos(2*PI(hx+ky+lz))
		  c = CosX[h]*CosKYplusLZ - SinX[h]*SinKYplusLZ;

		  s *= FiniteFactor[m][n];  // form factor * sin(2*PI(hx+ky+lz))
		  c *= FiniteFactor[m][n];  // form factor * cos(2*PI(hx+ky+lz))

      // add to real s.f. part
		  ReAmplitude[n] += c;  
      // add to imaginary s.f. part
		  ImAmplitude[n] += s;

      // store individual atom s.f. contributions to be used 
      // in TransformDerivatives()
		  v1[n][m]  = s;
		  v2[n][m]  = c;
		}
	}
}


